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Abstract 

Implementation of state variable-based viscoplasticity models is 
made in a general purpose finite element code for structural appli- 
cations of metals deformed at elevated temperature. Two consti- 
tutive models, i.e. Walker’s and Robinson’s models, are studied in 
conjunction with two implicit integration methods: the trapezoidal 
rule with Newton-Raphson iterations and an asymptotic integra- 
tion algorithm. A comparison is made between the two integration 
methods, and the latter method appears to be computationally more 
appealing in terms of numerical accuracy and CPU time. However, 
in order to make the asymptotic algorithm robust, it is necessary 
to include a self adaptive scheme with subincremental step control 
and error checking of the Jacobian matrix at the integration points. 
Three examples are given to illustrate the numerical aspects of the 
integration methods tested. 
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1 Introduction 


In the design of structural components, such as jet engines or nuclear pri- 
mary vessels, inelastic deformations and ratchetting effects under thermo- 
mechanical cyclings are of major concern in predicting the service lifes of 
such structures. To reliably predict the mechanical response of metals or 
composite materials at elevated temperatures, it appears that constitutive 
relations in the framework of a unified viscoplasticity theory are most ef- 
fective. Indeed, such constitutive models can represent typical deformation 
phenomena, e.g. creep-plasticity interactions, cyclic strain softening or hard- 
ening, and thermal ratchetting, none of which cam be properly described by 
the classical creep and plasticity theories. Although the viscoplastic theories 
are very useful in predicting the deformation behavior of high-temperature 
materials, they present considerable difficulties when the associated differen- 
tial equations are to be solved by numerical means. This is directly related 
to the ’’stiff" nature of the nonlinear differential equations, for which a small 
change of a state variable may result in large changes of others. 

In the literature, a great deal of study on numerical methods for integrat- 
ing viscoplastic rate equations in relation to the solution of boundary value 
problems has been reported [1-3]. Several numerical schemes have proven 
to be quite effective [4]. In general terms, there are one-step methods verse 
multi-step methods, such as those due to Gear [5]. It has been argued that 
although the multi-step methods generally give better solution accuracy as 
well as reliability, they are not suitable for large-scale finite element analysis 
because of their excessive demand in computer storage. Among the one-step 
integration schemes, two classes of numerical methods can be identified: ex- 
plicit and implicit schemes. Considerable discussion on the advantages and 
disadvantages of these two methods is available, e.g. [1,2]. No attempt will 
be made in the present paper to duplicate what is already available in the 
literature. Nevertheless, two observations may be offered: 1) the use of an 
explicit scheme in conjunction with self-adaptive time stepping can be very 
effective, however the solution accuracy may vary with the accuracy of error 
check embeded in the method, and 2) an implicit scheme together with iter- 
ations is generally more reliable, but requires intensive computational effort. 

More recently, an asymptotic integration algorithm was proposed by 
Walker, et al. [6-8], where the constitutive rate equations are cast into an in- 
tegral form. The integrands of solution are then expanded into Taylor series 
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in a time interval, say t G [t n ,t n + At], so that the resulting equations can 
be integrated term by term. When the series expansion is done about the 
upper limit [f„ + At], it gives an implicit iterative scheme and the numerical 
solution is unconditionally stable provided that the corresponding solution 
is exponentially decaying [7]. One apparent advantage of this method is 
that the numerical scheme involves the solution of only a 3 x 3 matrix of 
equations as opposed to a 13 x 13 (for 2D problems) or a 19 x 19 (for 3D 
problems) matrix of equations in the case of other implicit schemes, such as 
the trapezoidal rule with Newton-Raphson iterations [3,9]. From a computa- 
tional standpoint, the asymptotic integra tion algo rithm appears to be quite 
appealing. Nevertheless, the method is relatively new and it has not been 
tested on various viscoplastic models. One of the objectives of this paper 
is to conduct a comparative study between two, one-step, implicit methods, 
i.e. the asymptotic integration method and the trapezoidal method with 
Newton-Raphson iterations. The numerical schemes will be tested on two 
constitutive theories: Walker’s [10] and Robinson’s models [11] for isotropic 
materials. Moreover, another objective of the present paper is to show how 
the viscoplastic models and the associated integration schemes can be conve- 
niently implemented into a finite element code for structural analysis. Three 
numerical examples are included in this paper. 


2 Constitutive Models 

The present paper is concerned with the thermo-mechanical behavior of poly- 
crystalline metals and similar alloys deformed at elevated temperatures. In 
this connection, inelastic phenomena such as creep-plasticity interactions, 
cyclic strain-hardening or softening, and thermal ratchetting are of primary 
interest in a structured anedysis. Although a number of viscoplasticity theo- 
ries are available to represent the aforementioned phenomena, see e.g. [12], 
our attention is focused on a class of state- variable based viscoplastic models, 
that are based on a phenomenological approach, as discussed by Freed and 
Chaboche [13]. 

In the present study, the material is assumed to be initially homogeneous 
and isotropic. Furthermore, the analysis is limited to small deformations. 
Referring to a rectangular Cartesian coordinate system z,-, i=l,2,3, the total 
strain rate at a material point in a deformable body is given by 
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( 1 ) 


i = i e + e T +&T 

with the incompressibility condition in inelastic strains 


trU 1 ) = 0 (2) 

The relation between the elastic strain and the Cauchy stress rates is 
given by the Hooke’s law: 

d = D:(£-£ 7 -ar) (3) 

where £ is the total strain, £* is the elastic part of the strain, g is the matrix 
of thermal coefficients, T is the temperature , a is the Cauchy stress and D 
is the elastic stiffness matrix. 

To represent the inelastic process of a viscoplastic material, two sets of 
laws are being used: 1) a flow law which governs the rate of inelastic strain 
as a function of the current deformation state of a body, and 2) a set of 
evolutional equations that defines the rate of change in the internal state 
variables. The framework of these laws are outlined as below. 
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Flow Lslw : 


s' = /( §)S (4) 

where / represents a plastic scalar function, S is the effective deviatoric 
stress and a is the back stress. 


Evolution Equations : 

Evolution equations describe the change of internal state during the de- 
formation process. In general, three types of state variables can be iden- 
tified: back stress, drag strength and yield strength [9 ,10], At present, we 
shall consider viscoplastic models that do not have an explicit yield surface. 
Moreover, the change of state variables is described by two c omp eting mech- 
anisms: hardening vs. recovery processes (including both dynamic and static 
recovery). With the aforementioned considerations in mind, evolution equa- 
tions may be expressed for the back stress and drag strength, respectively, 
in the form: 


a = Hai 1 - R a a (5) 

k = H k \\i I \\ 2 -R k k (6) 

where a designates the back stress, a tensor, k is the drag strength, a scalr 
quantity, H represents a hardening function and R is a recovery function. 

For later discussion, both (5) and (6) may be combined into one repre- 
sentative form: 


= Hv.e 1 -Rn:Q (7) 

where Q = (a,fc). 

In addition, we define the following quantities: the deviatoric components 
of stress are 

§ = s- ^ tr (?) ( 8 ) 

The second invariants associated with the effective and the back stresses, 
respectively, are: 
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(9) 


h = : 5 ) 

h = ( 10 ) 

The material is assumed to behave thermo-elastically in hydrostatic response, 

tr(o-) = (A + | fi)[tr{i ) - 3aT] (11) 

where A is the Lame constant, fi is the shear modulus and a is the coefficient 
of thermal expansion. 

For computational purposes, we have selected two specific viscoplastic 
models for finite element implementation. These are Walker’s [9] and Robin- 
son’s models [10] for isotropic materials. The basic difference between the two 
models is their mathematical form. Walker’s model prescibes fairly smooth 
behavior in the state space. On the other hand, Robinson’s equations contain 
mathematical discontinuities in various regions of state space which require 
special care in numerical integration. Both models are briefly summarized 
below. 


2.1 Walker’s Model 

Flow law : 


i r w = f? 


with 


Evolution equations: 


f m T 


sfh 


k L k J 


n — 1 


( 12 ) 

(13) 


a(t) = (n a + n 2 )i I - a(t) - n x g 7 (t) (n 3 + n 4 e n6 )R + n e (-I 2 ) 2 




(14) 


k{t) = n 7 R - n 6 R k(t) - ko -n 9 [k(<)-ko] P (15) 
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In the above, R is the magnitude of inelastic strain rate defined by: 



and the effective stress is 

5 = fS-« (17) 

where n, n 2 , 113 , n 4, n s , n 6 ,nr, p, fco are material constants, which are 
defined in [10]. 


2.2 Robinson’s Model 

In order to capture very different viscoplastic responses of a material under 
different loading paths, several discontinuous functions were introduced in 
Robinson’s model, which was derived ftom an assumed plastic potential. 
This model constitutes an inelastic flow law and one evolution equation as 
follows. 

V 

Flow lag : 


where 


P(g)Mr) 

2rkT^/Ji 


(18) 


S = 5-a (19) 

In the above, P is a spline function defined in the stress or strain space [11], 
which provides a smooth transition from one set of material equations to 
another; F and G are material functions, also given in [11], 


Evolution equation: 


H a .j Ra$G{T) ^(m-0) 

a = ; t =— ( j v 'a 

~ G&~ k T y/T 2 


( 20 ) 
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where ifc r ,n,r,m,/?, R a and H a are material constants to be determined ex- 
perimentally and $(T) is a function of the temperature to account for tem- 
perature dependence of material properties [11]. 


3 Numerical Procedures 


With the constitutive equations outlined in the preceeding section, one may 
proceed to the structural analysis for viscoplastic materials using the finite 
element method. In this context, two solution approaches are possible: the 
initial strain method and the tangent stiffness method. In the former ap- 
proach, the viscoplastic deformations are treated as the initial strains on the 
right hand side of the incremental equilibrium equations. The major advan- 
tage of this method is that the global stiffness matrix remains constant under 
isothermal condition but the convergence rate of the corresponding iterative 
solution is at best linear. As an alternative, it is possible to derive the tan- 
gent material stiffness for a viscoplastic model [12,14,15] so that a tangent 
stiffness method can be pursued. In this context, the global stiffness matrix 
has to be reformed for each loading increment, but the iterative solution is 
quadratic in terms of its convergence rate. The, choice between these two 
methods is largely hinged upon the coding convenience. In our analysis, we 
have chosen the initial strain method. 

The finite element viscoplastic analysis generally involves two levels of 
numerical computations: global equilibrium and local constitutive calcula- 
tions. At the global (or structural) level, the targeted solution time is divided 
into a number of time intervals and the incremental equilibrium equations 
are solved for each time step in succession. At the local level (or material 
point), the constitutive rate equations are integrated by a numerical scheme 
where the stress and strain rates are converted into incremental quantities. 
As a result, the rate form of the constitutive equations in (3) for a typical 
time step t' (E [t,t + At] is written as 

Ac = D : (Ae - Ae 1 - aAT) (21) 

where the vector of incremental stresses are obtained from 



t+At 


adt 


( 22 ) 
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The expressions of Ae , At 1 and AT follow similar definitions. 

The main objective of our analysis is to integrate the rate equations (1) - 
(4) and (7) in an efficient way. Due to the ’’stiff” nature and the nonlinearity 
of viscoplastic rate equations, the solution accuracy of Ae 1 is very sensitive 
to the values of the state variables evaluated at the material sampling points. 
Therefore, an accurate and reliable numerical scheme is essential to integrate 
the aforementioned rate equations. 

For the sake of brevity, we write (l)-(4) and (7) in the condensed form 

?=/(*.'> < 23 > 
The above equation represents a system of nonlinear, first-order, ordinary 
differential equations, where y — (o-,a,£ J ,fc). For a three-dimensional prob- 
lem, y represents a vector of 19-components, i.e. 6-Cauchy stresses, 6 back 
stresses, 6 inelastic strains and 1 drag strength. For a 2D problem, y contains 
13 components. 

A number of numerical integration schemes are available to integrate (23). 
The choice of a particular scheme is dominated by three considerations: 1) 
suitability for finite element implementation, 2) solution accuracy and relia- 
bility and 3) computational cost. For instance, Gear’s multi-step methods [5] 
are known to be very effective for integrating nonlinear differential equations 
such as (23) for viscoplastic materials under homogeneous deformation. How- 
ever, these methods are not particular suitable for large-scale finite element 
computations for two reasons: 1) the methods require extensive central mem- 
ory and 2) a special start-up procedure is needed. In view of this, one-step 
integration methods are much more desirable. 

In the context of one-step integration methods, two classes of algorithms 
can be identified: explicit and implicit algorithms. The explicit algorithms 
(although less computation effort is required) are not very reliable, especially 
in dealing with thermo-mechanical cyclic loadings. For added solution ac- 
curacy, we have chosen the implicit algorithms. In particular, we select the 
trapezoidal rule with Newton-Raphson iterations and an asymptotic algo- 
rithm proposed by Walker [6,7]. Both algorithms are briefly outlined in the 
following subsections. 
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3.1 Implicit Trapezoidal Rule 

We consider a typical time interval re [t,t + At]. The trapezoidal rule can 
be easily obtained by using Taylor series expansion of the rates y between 
succesive iterations at time points t and t + At. For discussion purposes, we 
introduce the following notations: y n = ^(t),^n+i = Jf(^ "h A 0 correspond- 
ing to the i-th iteration. Thus, the trapezoidal rule with Newton-Raphson 
iterations prescribes [3,9,18,19]: 



r At^h 
r 2 dy J 

. A 

-y n y\+ 1+ 2 [/" +/n+i] 

(24) 

with the recursive relationship 





ySf. 1 ’ 

= vU\ + A / 

(25) 


where (df /dy) 1 is a Jacobian matrix evaluated at the time t-t- At for the i-th 
iteration The size of the Jacobian matrix is: 13 x 13 for 2D problems, and 
19 x 19 for 3D problems. 

It is known that stability of this algorithm is assured by its implicit- 
ness[3,9]. However, some comments are needed on the uniqueness of the 
solution scheme. In order to assure the uniqueness of Ay in (24), it is nec- 
essary to impose a condition that the p th norm of the Jacobian matrix must 
always be less than one. p may be any norm as long as it satisfies the basic 
matrix relationships [3,5,18,19]. 

We note again that (24) represents a system of unsymmetric simultaneous 
equations which have to be solved at every integration point of each finite 
element. For a sizable finite element mesh (say 1000 elements) in three- 
dimensions, the computation efforts required to solve these equations can be 
quite intensive. 

3.2 Uniformly Valid Asymptotic Algorithm 

The basic idea behind the algorithm proposed by Walker [6,7] is to trans- 
form the differential equations into a set of integral equations, which can 
then be solved approximately using a recursive relationship. In order to 
evaluate the resulting integral expression, an asymptotic expansion of the 
related integrand is performed about the upper limit of the time interval 
[t,t + At], resulting in an implicit integration scheme. The main advantage 
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of this method is that only a 2 x 2 (or for certain viscoplastic models, 3x3) 
matrix of equations need to be solved during the iteration process. 

After some manipulations, the viscoplastic rate equations, i.e. (l)-(4) and 
(7), can be written in the following symbolic form [6,7] 

y + X:y = H :g (26) 

^ * w 

where X and H are diagonal matrices and each contains only 2 or 3 distinct 
entries. The matrix H may be a function of time [10]. Within a typical time 
step [t,t + At], one can obtain the analytical solution of the above equation, 
i.e. 


y(t + At) = y(()e- iX + ■/■<*" H(()j({)<« ( 27 ) 

~ ~ J(—t ** 

where 


AX = X(< + At) - X(t) « X • At (28) 

Here y(t ) is the initial value of the appropriate stresses and strains defined 
at the begining of the current step. The function g is defined for a specific 
material and is, in general, a function of the current state and the appropriate 
rates. 

We denote the integral expression on the right hand side of (27) by I (At) 
so that 


y(t + At) = e -ax ■ y(t) + I(At) (29) 

and 

/(At) = jT 4 ' e-W‘ +il »- x «»H(o|({K (30) 

It is seen from the above, in order to solve for y(t + At) the integral I (At) 
has to be evaluated, and the entries in X and in g still remain unknown. 
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Approximation of Integral : 

If the solution is exponentially decaying, it is possible to approximate the 
integrand in (30) by an infinite series expression 
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I (At) = jf e~ ax JJ J 4„z n ]<iz (31) 

Thus, by retaining the first two powers of ”z”, we obtain the following ap- 
proximate expression for the i ^ component of the integral 


-X:At 


Ii(At) = Hi&Agi 


1 - <>-•*•** . e _x * At - ■ j~ rr 

+ «*<)’ ^ 

XiAt 


XiAt 


,no sum over 1 

(32) 


where the derivatives of g are approximated using the known values of this 
vector at the begining and at the end of the current step: 


gi(t + At) - gj(t) 
At 


(33) 


Equation (32) is similar to the one derived in [6-8]. Here it incorporates the 

time dependency of the matrix H [10]. 

With the numerical equations outlined in the general form, it is useful to 
list the specific relations for the physical quantities, i.e. back stress, effective 
stress and inelastic strain. 

Back stress: 


Effective stress: 


Inelastic strain: 


a(t + At) = exp AX a(f) -f I a (At) 

(34) 

E(t + At) = exp -AX S(t) + 7 s (At) 

(35) 

: , = e(( + A() _5(i±^)±*(i±M 

(36) 


where e is the deviatoric strain tensor. 


Evaluation of Matrix AX : 
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The entries of diagonal matrix AX appeared in (26) are related to the 
second invariants of inelastic deformations, which are generally unknown. 
These quantities are determined by an iterative process outlined as follows. 
First of all, we recognize that the matrix AX must be semipositive, i.e. 

AX > 0. (37) 

Since AX is a diagonal matrix, for convenience it can be treated as an equiv- 
alent vector in its diagonal form. From (7), the general form of AX may be 
written as 


AX = R [y(t + A<)]A< (38) 

This vector equation defines an implicit nonlinear system which can be 
solved by an iterative procedure. To this end, we define a residual vectorial 
function F 


F = AX - R(y)Ai (39) 

mg - ... ........ ..... 

The above equation is solved by a Newton-Raphson iterative procedure for 
the incremental change of AX, which is denoted by A AX. Invoking Taylor 
series expansion and ignoring higher order terms, we find 


i-( 


dR 

SAX 


)« 


AAX = -F (i > 


(40) 


where AAX represents an error vector between two successive iterations. 
The updated value for AX is thus given by 

AX (i+1 > = AX (i) -I- AAX (41) 

Iterations are terminated when the following convergence criterion is satisfied 


11 A AX|[ 2 
I|ax|| 2 - c 


(42) 


where £ is a tolerance limit and || • || 2 designates a Euclidean norm. It is 
noted again that (40) represents a system of 2 x 2 simultaneous equations re- 
gardless whether for a two-dimensional or three-dimensional problems. This 
is obviously a major advantage of the asymptotic integration algorithm over 
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the trapezoidal rule where one has to solve a system of 19 x 19 simultaneous 
equations for every iteration. 

In (40), the coeffficient matrix, denoted by J = represents a 2 x 2 
Jacobian matrix, which may be derived analytically or defined by a numer- 
ical means. To improve computational efficiency of the procedure, we have 
obtained the analytical expressions for both Walker’s and Robinson’s models. 

3.3 Time Step Control 

As mentioned in the previous susections, numerical stability of trapezoidal 
and asymptotic integration algorithms is assured by the fact that both pro- 
cedures are implicit methods. However, special care must be given to insure 
solution uniqueness, or the invertibility of equation (40). To this end, a 
Lipschitz condition [3,5,18] is employed. The Lipschitz number is defined by: 

L = sup ||J(OIIp ( 43 ) 

£ eAi 

where || || p is the p th order norm. 

Lipschitz theorem states that if 0 < X <. 1 then there exists a unique 
solution and the iterates AX’ will converge after sufficient number of ita- 
rations is performed. This theorem defines a convergence condition for the 
Jacobian matrix, i.e. 


max |Ai(J)| < 1 1 J | |p < 1 (44) 

i 

where A, is the i-th eigenvalue of J. The requirement in (44) insures the 
existence of the inverse of (I — J). The particular p-norm is at the discretion 
of the user. Fulfilling this condition, the algorithm is unconditionally stable 
and will yield a unique solution. 

Our time step control makes use of the Lipschitz number. That is, dur- 
ing the analysis the magnitude of the Jacobian matrix is being monitored. 
Whenever its norm exceedes an allowable limit, the global time step is subdi- 
vided by the calculated Lipschitz number. Thus, the subincrement step size 
is given by 


At = 


At 

m 


(45) 
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where L(t) is the current Lipschitz number and At is the global time step. 

4 Computer Implementation 

The viscoplastic models discussed in the previous sections have been written 
in the form of a material module for finite element computations. Although 
the module was written for implementation into a specific code, i.e. NFAP 
[16], it is adaptable to other general purpose finite element codes with mini- 
mum effort. . r , _ _ . . •... .. -- i! ^ l_. 

Basically, the material module constitutes two main branches: integra- 
tion methods and material functions. In the integration branch, one has 
the option to choose explicit schemes, e.g. forward Euler or Runge-Kutta 
methods, or implicit schemes, such as the trapezoidal rule or the asymptotic 
integration method. In the material function branch, one may select Walker 
[10], Robinson [11] or any other viscoplastic model. The coding is transpar- 
ent enough so that user may easily add another integration method and/or 
viscoplastic model. 

The material data are calculated at all integration (or Gauss) points of 
all elements: that is, all continuum elements (2D, 3D, plate or shell ele- 
ment) share the same material module. The material module is called by the 
main program in three stages for each global time step: 1) element stiffness 
calculation, 2) calculation of the internal forces for global equilibrium check 
(iterations), and 3) stress and strain output. Data transfer between the main 
program and the material module is achieved through the subroutine argu- 
ment list and common blocks. In entering the material module, the following 
data are defined: £ n > |Jn) T* n , fnj At^ wid where ( ) n 

indicates quantity at time t; ( ) n+1 , quantity at time [t-f At]. When calcula- 
tions are done in the material subroutines, the material stiffness matrix and 
the updated inelastic strains, stresses and state variables corresponding to 
the current deformation are passed onto the main program. A skeletal block 
diagram for the material module outlining the major calculation controls are 
shown in fig.(l) and specific functions of each control routine are described 
below. 


1. Main driver - Updates the solution vector y , forms the material stiffness 
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matrix D(f + At), and if necessary, adjusts the global time increment 
At to obtain the local subincrement At at each Gauss point. 

The main driver calls upon two control routines: numerical algorithm 
control and subincremental time control. 

2. Numerical algorithm control - At each time step, it selects the appropri- 
ate algorithm optioned by the analyst, and computes the corresponding 
iteration vector y' n+ i with a suitable convergence test. 

3. Material function control - Calculates the inelastic strain and associated 
evolution equations for a specific viscoplastic model. If the asymptotic 
algorithm is employed, it evaluates the coefficient matrices X, H and 
the vector g defined in (26) for either Walker’s or Robinson’s model. 

4. Jacobian matrix control - For implicit algorithms, it evaluates the Ja- 
cabian matrix of the material model and checks a criterion for solution 
convergence and invertibility of the matrix equations involved. 

5. Subincremental time control - Using the criterion defined in the previ- 
ous section, it determines whether or not the global time step A t is to 
be divided into subincrements: At = At/ L, where L is the Lipschitz 
number, see (45). 

5 Numerical Examples 

Three sample problems are included to demonstrate the utility of numerical 
algorithms implemented for finite element applications of viscoplastic struc- 
tures. Although the constitutive equations incorporated in the program can 
be utilized for any general three-dimensional state of stress, the problems 
considered here are merely two-dimensional. In our analysis, some numerical 
difficulty was encountered due to the mathematical discontinuity existed in 
Robinson’s constitutive equations and ways to alleviate the aforementioned 
difficulty are discussed. Further, a comparison of the numerical performance, 
in terms of CPU time and solution accuracy, between the trapezoidal rule 
and aymptotic algorithm is also included in the present study. 
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5.1 A Thick Walled Cylinder 

The first example is a thick walled cylinder subjected to a time dependent 
internal pressure, p(t), which undergoes a full loading- unloading cycle. The 
finite element model contains 12 axi-symmetric elements and is shown in 
fig.(2) along with the input pressure function. The material for the cylinder 
is represented by Walker’s isotropic model, which does not have any diconti- 
nuities in its mathematical formulation. The main objective of solving this 
problem is to verify the coding for the various numerical integration proce- 
dures considered. - 

Material constants for this problem are [6]: 

T = 1600 F E = 18.6 • 10 3 Ksi ko = 91 505 Ksi v = 0.345 

p =3. n =5.128 n, = 1.158 n 2 = 0. 

n 3 = 5 • 10® n 4 = 672.6 n B = 0 n 6 = 0 

n 7 = 8.98 • 10 -4 n 8 = 0 n 9 = 0 

To carry out the incremental analysis, the entire load cycle is divided into 
50 time steps with variable step sizes within the four quarters of the load cy- 
cle. Solution was obtained by using four different integration schemes: the 
asymptotic algorithm, the trapezoidal rule with Newton-Raphson iterations, 
the forward Euler and Runge-Kutta procedures with error checks [3,9,12,18]. 
The resulting hoop stresses are plotted in figs. (3) and (4). As seen in the 
plots, no visible difference can be detected among the numerical results ob- 
tained by the different integration schemes. 

The CPU time for each of the integration schemes is listed in Table 1. It 
appears that both the Runge-Kutta method and th e asymptotic algorithm 
consume about 20% less CPU time than the implicit schemes. The Euler 
method seems to be the most time-consuming algorithm in order to achieve 
the same degree of solution accuracy. 


5.2 A Simply Supported Beam 

The second example is a simply-supported beam under a state of plane stress. 
The loading function is a full-cycle displacement prescribed at the center of 
the beam. The finite element model and the applied load function are plot- 
ted in fig. (5). For this problem, Robinson’s isotropic model is employed. As 
noted before, Robinson’s model contains mathematical discontinuities which 
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may pose numerical diffculty during the course of integrating the constitu- 
tive equations. We use this problem to illustrate how numerical difficulty 
may arise as a result of discontinuities existing in the constitutive equations. 
To remove such difficulty, we need to monitor the variation of the Jacobian 
matrix so that its norm must be kept within a certain allowable tolerance. 
Otherwise, numerical anomaly will occur, which in turn causes either erro- 
neous or non-unique solutions. 

The material constants used for this example are[10]: 


T = 1000.4 F 
n =4 
R = 9 • 10- 8 
W x = 0.1 


E = 22,480 Ksi 
r = 3.61 • 10 7 
H = 1.37 -10 4 
W 2 = 0.1 


k T = 0.82 Ksi 
m = 7.73 
G 0 = 0.1 


v =0.345 
P = 1.5 


Regions of Numerical Difficulties - When solving the beam problem using 
Robinson’s model with the asymptotic algorithm, oscillations in the numeri- 
cal solution were detected. There are in fact two different types of numerical 
oscillations: case 1) ||J|| > 1, and case 2) ||J|| ^ 0. The first case occurs 
when the rate of change of the state variables is very rapid (or high rates are 
encountered), thus the corresponding constitutive rate equations are in a stiff 
regime. The second case represents another extreme; that is, the material 
is almost responding elastically with very little inelastic deformation. Thus, 
solution of (26) leads to numerical drifting. 

To further illustrate the numerical phenomenon stated above, we have 
plotted in fig.(6) the bending stress vs. the back stress at an integration point 
near the upper fiber near the center of the beam. As seen in the figure, the 
back stress is changing much faster then the corresponding stress component, 
thus leading to singularity of the coefficient matrix in (26). This phenomenon 
is further verified from the plots of the Jacobian norm vs. the number of 
local iterations at the integration point in fig.(7). The end result is that a 
considerable amount of CPU time is consumed for very little improvement 
in solution accuracy. In fig. (8), the dots close to the abscissa represent the 
drifting of the Jacobian matrix. To alleviate such difficulty, we take an 
average value of AX in (40) evaluated between two successive iterations. By 
doing this, the numerical drift was immediately elimianated. 

To stabilize the numerical oscillation of the first case, a step control pro- 
cedure was implemented that checks the magnitude of Jacobian norm. When 
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the norm of the coefficient matrix in (40), ||I — J|| is kept within a tolerance 
limit, oscilations in the numerical solution of (26) are dramatically reduced 
to relatively stable values as indicated in figs.(9). It shows a reduction in 
the Jacobian norm when the integration steps are limited within a reason- 
able value. As a result, the integration algorithm yields a unique convergent 
solution. 

The same problem was also analyzed for the trapezoidal rule and Runge- 
Kutta method with error control [17]. In the case of trapezoidal rule, a 
constant fixed-step history was employed, whereas an adaptive step proce- 
dure was activated in the case of asymptotic algorithm. The trapezoidal rule 
failed to yield a convergent solution in the critical discontinuity region as 
indicated in fig. (10). On the other hand, the solutions obtained from the 
asymptotic algorithm and the Runge-Kutta method converged quite nicely, 
since both methods have a built-in step control. 

The CPU time for the three integration methods is detailed in Table (2). 
The time required by the asymptotic algorithm is about 50% less than that 
of the trapezoidal rule. For a larger size problem, the time saving will be 
much more significant. 

5.3 A Cylindrical Thrust Chamber 

The third example is a cross section of a cylindrical rocket thrust chamber 
which is subjected to thermal as well as pressure cycling [19]. Due to the 
repetitous pattern of its cross section, only one sector including the cooling 
channel is considered for finite element analysis, fig. (11). In addition to 
the thermal loading, the chamber is subjected to a time dependent internal 
pressure; p\(t) acting on the cooling channel surface, and p 2 {t) acting on 
its inner surface. The temperature and pressure functions p x and p 2 follow 
the same time history and the corresponding function for one typical loading 
cycle is shown in fig. (12). The sector of the chamber is modeled by 34 8-node 
plane strain elements, fig. (12). The chamber material is assumed to follow 
Robinson’s isotropic model with the same material constants as were used 
for the beam. In the finite element analysis, both the asymptotic algorithm 
and trapezoidal rule were used for comparison of solution efficiency. 

Shown in fig.(13)is the deformed shape of the thrust chamber cross section 
at the end of one loading cycle. From the analysis, considerable inelastic 
strain accumulations occurred at the inner corner of the cooling channel, 
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which in turn caused bowing effect near the central section. This is called 
the ”dog house” effect and it is observed in experiments [20]. Also shown in 
fig.(14) is the maximum inelastic strain history that occurred at the corner 
of the cooling channel. 

We use this problem to compare the solution efficiency between the asymp- 
totic algorithm and trapezoidal rule. To this end, we employ two different 
step histories: 1) a constant step history for both methods as a basis of com- 
parison, and 2) a variable step history with step control for the asymptotic 
algorithm. After several trials, it was found that 1000 constant, global, time 
step increments were needed to obtain a convergent solution for the trape- 
zoidal integration rule. Table (3) summarizes the CPU time ratios between 
the two methods. The CPU time ratio per iteration is defined as 


C PU Ratio = 


CPU time using trapezoidal rule 
CPU time using asymptotic algorithm 


As seen in the table, the ratios vary somewhat during the course of the 
analysis, i.e. within the range of 1 — 5. It is pointed out that the asymptotic 
algorithm is particularly efficient when the rate of change of the state vari- 
ables becomes large. In addition to the constant ( step history, we employed 
the asymptotic algorithm with an adaptive step control. In this case, the 
method becomes even more attractive since it requires only about one-fourth 
of the CPU time consumed by the trapezoidal method. 


6 Conclusion 

Two state variable-based viscoplastic models, namely Walker s and Robin- 
son’s models, have been implemented into a general purpose finite element 
code for structural applications of metals deformed at elevated temperatures. 
These models are represented in a material module form so that they may be 
easily modified or transplanted to other finite element codes. The addition 
of new viscoplastic models or integration routines is relatively transparent 
to an investigator. A comparative study has been made on the numerical 
performance of two implicit integration methods: the trapezoidal rule and 
the asymptotic algorithm. It was found that the asymptotic algorithm be- 
comes very effective if a self adaptive scheme based on an error check of the 
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Jacobian matrix and subincremental step control is utilized at the material’s 
integration points. 
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Table 1: CPU Ratio of a Thick Walled Cylinder 
( Walker’s Model ) 



Explicit 
Runge Kutta 

Implicit 

trap. 

Forward 

Euler 

Cal. of 
Nonlin. 
Mat. 

1.168 

1.137 

1.7 


Table 2: CPU Ratio of a Simply Supported Beam 

( Robinson’s Model ) 



Explicit 
Runge-Kutta 
25 steps 

Explicit 
Runge-Kutta 
50 steps 

Implicit 

Trap. 

Cal. of 
Nonlin. 
Mat. 

1.8 

2.5 

no 

convergence 


Table 3: CPU Ratio of a'Thpist Chamber 
( Robinson’s Model ) 


Step 

Range 

CPU 

Ratio 

1- 2 

5 

3- 200 

1 

201- 400 

2 

401- 600 

2.5 

601- 800 

2 

801-1000 

3 
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Figure 1: Block Diagram of Material Module 
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Figure 10: Robinson’s Isotropic Material for a Simply Supported Beam 







Figure 11: Cross Section of a Cylindrical Thrust Chamber 
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Figure 13: Deformed Geometry of the Cylindrical Thrust Chamber 
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Figure 14: Strain History (e-*) at Point A in the Thrust Chamber 
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